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ABSTRACT 

The contour integral method of Hunter & Qian is applied to axisymmetric galaxy 
models in which the distribution function (DF) is of the form / = f{E, L^), where E 
and Lz are the classical integrals of motion in an axisymmetric potential. A practical 
way to construct the unique even part fe{E,Lz) of the two-integral DF for such 
On systems is presented. It is applied to models, both oblate and prolate, in which the 

^ mass density is stratified on similar concentric spheroids. 

^ The spheroids with scale-free densities are discussed in detail. These provide 

useful approximations to the behaviour of more realistic models in the limit of small 
and large radii. The self-consistent case is treated, as well as the case in which there 
are additional contributions to the potential from a central black hole or dark halo. 
" ' The two-integral DFs for scale-free densities in a Kepler potential are particularly 

simple. These can be used to model power-law density cusps near a central black 
r hole, or to model the outer parts of finite-mass systems. The range of axis ratios 

^*>0 and density profile slopes is determined for which spheroidal power-law cusps with a 

central black hole have a physical two-integral DF. 

More generally, the two-integral DFs are discussed for a set of spheroidal '(a,/3)- 
models', characterized by a power-law density cusp with slope a at small radii, and 
a power-law density fall-off with slope a + 2/3 at large radii. As an application, the 
DF is constructed for the (a,/3) model with a 1.8 x 10® black hole used by van 
^ der Marel et al. to interpret their high spatial resolution spectroscopic data for M32. 

The line-of-sight velocity profiles are calculated. The results confirm that the model 
^ fits the data remarkably well. The model is used to calculate the kinematic signatures 

?-H of a central black hole in observations such as are now possible with the Hubble 

(/3 Space Telescope. The predicted Gaussian velocity dispersion for the M32 centre is 

C3 127 km s-i with the 0.09" X 0.09" square aperture of the Faint Object Spectrograph, 

and 105 km s"-*^ with the 0.26" diameter circular aperture, while the central dispersion 
measured from ground-based data is only SGkms"-*^. 

Key words: stellar dynamics - galaxies: kinematics and dynamics - galaxies: struc- 
ture - galaxies: central black holes - galaxies: individual: M32 
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1 INTRODUCTION 

Few realistic dynamical models have been constructed for 
galaxies that are not spheres or discs. The main reason 
for this paucity of models is that in axisymmetric or triaxial 
galaxies the stellar motions are governed by non-classical in- 
tegrals of motion, which are generally not known explicitly. 
An exception is provided by axisymmetric models in which 
the phase space distribution function (DF) / = f{E, Lz), so 
that it depends only on the two classical integrals of motion, 
the energy E and the angular momentum component Lz 
parallel to the symmetry axis of the system. Hunter (1977) 



showed how the velocity dispersions in such a two-integral 
axisymmetric galaxy model can be calculated by solving the 
Jeans equations. Various authors have applied his solution 
to model kinematic observations of elliptical galaxies (e.g., 
Binney, Davies & lUingworth 1990; van der Marel, Binney & 
Davies 1990; van der Marel 1991; Cinzano & van der Marel 
1994; CaroUo & Danziger 1994). One disadvantage of this 
approach is that it is not evident whether the intrinsic veloc- 
ity dispersions that best fit the line-of-sight measurements 
indeed correspond to a physical model, in which / > 0. 

It is now possible to extract not only the mean line-of- 
sight velocity {■uios} and velocity dispersion irios from absorp- 
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tion line spectra, but also the shape of the line-of-sight ve- 
locity distribution, hereafter referred to as the velocity pro- 
file (VP) (e.g., Franx & lUingworth 1988; Bender 1990; Rix 
& White 1992; van der Marel & Franx 1993). To model 
such data one could continue to solve the Jeans equations of 
increasing order (e.g., Magorrian & Binney 1994), but it is 
preferable to calculate the entire DF, so that the theoreti- 
cal VPs can be calculated accurately, and only models with 
/ > are considered. 

The calculation of f{E, Lz) for axisymmetric models 
has long been hampered by certain perceived technical dif- 
ficulties, with as main result that only a few such DFs were 
found by various integral transform methods (e.g., Lynden- 
Bell 1962; Hunter 1975; Dejonghe 1986), usually for special 
mass models (but see Evans 1993, 1994; Evans & de Zeeuw 
1994). Hunter & Qian (1993, hereafter HQ) showed how 
these difficulties can be circumvented, and developed a con- 
tour integral method that in principle allows calculation of 
f{E, Lz) for a wide variety of mass models. In particular, 
it is no longer required that the density p(^B?,z^) can be 
written explicitly as a function of $ and B? , where $ is 
the relative gravitational potential (cf. Binney & Tremaine 
1987). 

In this paper we demonstrate how the HQ method 
can be used to calculate f{E, Lz) for realistic axisymmetric 
galaxy models, with emphasis on models stratified on simi- 
lar concentric spheroids with arbitrary density profiles. As a 
specific application, we study a set of '(a, /3)-models', char- 
acterized by a power-law density cusp with slope a at small 
radii, and a power-law density fall-off with slope a + 2/3 
at large radii. The HQ method allows us to include an ex- 
ternal potential, such as that of a central black hole or a 
dark halo. The present paper complements recent work by 
Dehnen & Gerhard (1994) who consider two-integral ax- 
isymmetric models with central density cusps, based on a 
convenient series expansion of f{E, Lz), and give a thor- 
ough discussion of the VPs. We derive the DFs for a larger 
class of models, include the effects of a central black hole or 
dark halo, and derive some further properties of the VPs. In 
addition, we present a detailed study of the spheroids with 
scale-free mass densities, and show how they can be used 
to approximate the dynamical structure of the more general 
models at small and large radii. 

We illustrate our technique by applying it to the nearby 
elliptical galaxy M32, which may contain a black hole (or at 
least a dark mass concentration) in its centre. Van der Marel 
et al. (1994b) used axisymmetric two-integral models to in- 
terpret the high spatial resolution kinematic observations of 
M32 by van der Marel et al. (1994a). The modelling con- 
sisted of: (i) use of Evans' (1994) power-law model DFs 
without a central black hole; and (ii) calculation of the first 
three moments of the VP for the case with a black hole, by 
solution of the moment equations of the coUisionless Boltz- 
mann equation. A remarkably good fit was obtained with a 
1.8 X 10^ Mq black hole, but the actual DF could not be cal- 
culated. With the technique presented here we can calculate 
the entire f{E, Lz) for the model with a central black hole, 
allowing a better comparison with the available data. We 
confirm and strengthen the results of van der Marel et al., 
and use the DF to calculate the kinematics and VP shapes 
that one would expect to observe with the high spatial res- 
olution of the Hubble Space Telescope (HST). 



This paper is organized as follows. In Section 2 we dis- 
cuss the implementation of the contour integral method to 
cases where the relation p = p(^^,B?) is known implicitly, 
show how for similar concentric spheroids f{E, Lz) can be 
found as a numerical quadrature for each E and Lz, and 
summarize how the VPs of two-integral axisymmetric mod- 
els can be calculated. In Section 3 we describe the properties 
of a family of spheroidal mass models with a central density 
cusp. We consider the special case of scale-free spheroidal 
models in detail, and discuss the inclusion of a central black 
hole. We apply the results to M32 in Section 4, and sum- 
marize our conclusions in Section 5. 

2 TWO-INTEGRAL DISTRIBUTION 
FUNCTIONS 

A general algorithm for the application of the HQ method 
to models in which the density as a function of potential 
and cylindrical radius is only known implicitly is presented 
in Section 2.1. The case of spheroidal mass models with 
arbitrary density profiles is discussed in Section 2.2. Section 
2.3 addresses the calculation of VPs for two-integral DFs. 
The reader who is interested mainly in the applications of 
the method can skip to Section 3. 

2.1 The HQ contour integral method 

We consider an axisymmetric model of infinite extent with 
a density p(^B?,z^) and an overall potential ^(^B? , z^). In 
a self-consistent system p and $ are related via Poisson's 
equation while in a non-self-consistent system $ contains 
contributions from other components, which may include a 
dark halo and/or a central black hole, besides that from the 
density p. The HQ method is applicable in both cases. We 
adopt the convention in which the potential attains its max- 
imum value at the centre and decreases outwards. Hence z^ 
is determined uniquely by the values of $ and R'^ , provided 
that the former lies between the potential at infinity and 
the equatorial potential $(il^,0). Therefore the density as 
a function of $ and R'^ , which we denote as p[^,R'^), can 
be obtained. This function, whose analytic continuation is 
needed in the HQ method, is implicit in cases where z^ can 
only be determined implicitly for a given pair ($, R'^). It is 
this implicit case of the contour integral method with which 
we shall be concerned. 

The HQ method can be used for both finite and infinite 
mass systems. When applied to the density p(^R'^ , z^) in a 
potential $(il^, z^) it gives the unique fe{E, Lz) that is even 
in Lz and that generates p. When applied to Rp{v^)[R'^ , z^) 
it gives the unique fo{E, Lz) that is oddin Lz and generates 
the mean azimuthal streaming motions (v^). In practice it 
is not easy to obtain good observational data on the two- 
dimensional mean {■uios } on the plane of the sky from which 
the intrinsic azimuthal mean streaming field {v^)[R, z) must 
be found. Therefore, an alternative approach is to take the 
odd part fo as a product of the even part /e and a prescribed 
function whose magnitude is no greater than unity. This 
ensures that / = /e + /o is physical (non-negative) if /e is. 

The physical values of E and Lz are those which cor- 
respond to bound orbits in the potential They lie in 
the region of the [E, £2)-plane that is bounded by the two 
straight lines Lz = and E = $<x>, and by the curve £, as 



Axisymmetric galaxy models with central black holes 3 



Figure 1. The physical domain of bound orbits in the 
(£/, L^)— plane (Lindblad diagram) is bounded by Lz = 0, 
E = ^oo and the locus £ of circular orbits in the equatorial 
plane, defined by equation (2.1). The thin straight lines through 
the point (£/,L^) indicated by the solid dot are tangent to £, and 
intersect the boundary Lz = at values ^min and ^max which 
bound the window V of physically achievable values of the po- 
tential energy ^ for orbits with integrals E and Lz- The special 
value ^env in V is the intersection of the straight line that is 
tangent to £ at energy E and Lz equal to the maximum allowed 
value Lc- 

shown in Figure 1. When = — oo, this region extends 
indefinitely to the left. The curve £ is the locus of the cir- 
cular orbits in the equatorial plane. It is described by the 
parametric relation 



E = ^iRl,o) + R; 



d$(J^^o) 



-2R 



dR^ 



(2.1) 



where Rc is the radius of the circular orbit. The value 
Lc = Lz{Rc) is the maximum allowed value of Lz at fixed 
energy. The set of orbits with energy E and angular momen- 
tum Lz cover a range [$min, ^max] of physically achievable 
values of the potential energy This range can be found 
geometrically by constructing the two straight lines through 
the point (^E, i^) that are tangent to the locus £ in the Lind- 
blad diagram (Figure 1). Their respective intersections with 
the boundary Lz = give $max and $min, and so mark a 
window V on this line (HQ). We denote the potential energy 
of the circular orbit of energy E by $env(-B). Its value can 
be found geometrically (Figure 1), and it also follows easily 
upon solution of the first of equations (2.1) for Rc, and then 
using $env(-B) = ^(-RciO). The window V and the value 
^env(-B) in it are important in the evaluation of f{^E, Lz) by 
means of the HQ method. 

One way of writing the HQ solution for the even part 
f^{E, Lz) of the DF is 



ME,Lz)- 



[*=„(B)+] 



47r2i^/2 



d^ 



(2.2) 



This is a complex contour integral on the complex ^-plane. 
The "density" term of the integrand is now a function of the 
single complex variable ^ (the complex potential), and the 
two subscripts 1 denote the second partial derivative with 
respect to its first argument (as in HQ). The value of this 
function in the complex domain is obtained via the analytic 
continuation of the physically relevant value /5ii($,il^). For 
simplicity we also denote the analytic continuation by pn. 
The physically achievable values of $ lie on the real ^ axis to 
the right of the point ^ = £^ in the window V, i.e., for ^ = $ 
in V, the point [^,R^ = |i^/(* - E)] lies in the physical 
domain of pn. Obviously values of pn on the window V 
must coincide with the physically relevant values there. For 
the square root term in the integrand we choose the branch 
induced by a cut to the left of the point ^ = E along the 
real ^-axis so that it is real and positive when ^ > E. This 
choice together with the fact that V is to the right of ^ = E 
ensures that the integrand is real for ^ in V. 

The point ^ = $env(-B), which only depends on E, is 
a point that always lies in the window V. As indicated by 
the notation [$env(-B)+] in equation (2.2), the path of the 
complex contour integral is taken as a loop which starts on 
the lower side of the real ^-axis at ^ = $<x>, crosses the 
real ^-axis at ^ = $env(-B), and ends at ^ = on the 

upper side of the real ^-axis. To evaluate the integral (2.2) 
we first specify a contour. Depending on whether is 
finite or (negative) infinite, a different parametrization of 
the contour must be used (Figure 2). For finite, one 
can always choose = 0. A simple parametrization which 
corresponds to an elliptical path is given by the equation 



-TV <e <TV. (2.3) 



-■K <e <-K. (2.4) 



2{$-E)\ ii-Efl^' 



i = |*env(-B)(l + cos6l) +i/isin6l. 
For = — oo we can choose 

Q 

^ = $env(-B) + 1{1 - sec -) + i /i sin 9, 

The parameter h determines the maximum width of the con- 
tour in both equations (2.3) and (2.4). It must be kept small 
in order for the contour to avoid enclosing any complex con- 
jugate singularities that the "density" term pn might have. 
However too small a value of h will force the contour to come 
near enclosed singularities on the real ^-axis, when greater 
care must be taken to achieve good accuracy of the numer- 
ical integration. The parameter I in equation (2.4) allows 
the location of the points of maximum width of the contour 
at 6 = ±7r/2 to be adjusted. We take both h and I to be 
of the order of 0.1$env(-B). These parametrizations convert 
the complex contour integral (2.2) into an integration with 
respect to the angle 6, and have been satisfactory in our 
computations. Other parametrizations are possible though. 
For = —00 it is sometimes convenient to first change 
the integration variable in the solution (2.2), to obtain an 
integral with a finite path (see Appendix B). 

The fact that the integrand in equation (2.2) is real- 
valued on V has a useful consequence. According to the 
Schwarz Reflection Principle (Levinson & Redheffer 1970), 
once we have succeeded in continuing the integrand into a 
domain above the real ^-axis, we can also continue it ana- 
lytically as a complex conjugate into the reflected domain 
below the real ^-axis. Therefore we can evaluate /e by inte- 
grating along either the upper or the lower half of the loop 
and multiplying the result by a factor of 2. For definiteness 
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we use the value at the previous point as the initial guess 
and employ Newton's method to solve equation (2.6) iter- 
atively. Numerical calculations we have done so far show 
that this approach provides close initial guesses for so 
that the iterations quickly converge to the correct branch. 
Once z^ is found we use equation (2.5) to evaluate the in- 
tegrand, which is then added to the weighted sum for the 
Gauss quadrature. As we approach ^ = along the con- 
tour, B? = ^L'l/(^^ — E) approaches a finite value, hence z^ 
becomes unbounded. Consequently any numerical method 
that we use for equation (2.6) breaks down. However the 
"density" term pn is vanishingly small at large distances 
for centrally condensed systems, so that its contribution to 
the integral (2.2) becomes negligible. This can be measured 
by the ratio of the value of the integrand to the value of 
the weighted sum at each given quadrature point. Once this 
ratio falls below a preset tolerance we stop and accept the 
weighted sum as the value of the integral (2.2). 



2.2 Classical spheroids 



Figure 2. The complex ^— plane, with the contour used in the nu- 
merical evaluation of f[E,Lz) by means of equation (2.2). Top 
panel: the contour (2.3), for the case when ^oo = 0. Bottom 
panel: the contour (2.4) for the case when ^oo = — oo. The win- 
dow V of physically allowed values of ^ is indicated by the thick 
solid bar along the real ^— axis (see Figure 1). The contour in- 
tersects this window at the value ^env (■£/), which is the potential 
energy of the circular orbit with energy E. 

we shall use the upper half for our calculations. 

To evaluate the integral (2.2) numerically for a given 
pair [E, Lz), we first discretize the contour (2.3) or (2.4) 
by a Gauss quadrature in 6, and then approximate the in- 
tegral by a weighted sum of the values of the integrand at 
the quadrature points. The main task is then to evaluate 
pii[^,B? = ^L'l/(^^ — E)] at each quadrature point ^ on the 
contour. By implicit differentiation, we obtain 



/5ii(^,-R') 



P22{R\z^) P2{R\Z^)^22{R\Z^) 



[^2{R^,Z^)]^ 



[^2{R^,Z^)Y 



(2.5) 



in which each subscript 2 denotes a partial differentiation 
with respect to z^ . In using this equation, we only have to 
find the value of z^ for a given pair = ^L'l/(^^ — E)]. 

For the case under consideration, a complex equation has to 
be solved at each quadrature point. We emphasize that the 
contour integral solution requires that the integrand attains 
its physically achieved values in the window V and its values 
on the complex contour are from the analytic continuation 
of the density. Therefore it is absolutely essential that the 
values of z^ , which are to be obtained from the numerical 
solution of the equation 

satisfy this requirement. 

We thus proceed as follows. We choose a pair [E, Lz), 
and start at the point ^ = $env(-B) {d = 0) in the window 
V. Since we are now in the physical domain, we look for the 
unique real positive solution of z^ of equation (2.6). As we 
move along the contour and arrive at a quadrature point. 



We now restrict our attention to axisymmetric systems of 
infinite extent in which the density of stars is stratified on 
similar concentric spheroids, i.e. 



p{m^), 



R^+z^/q\ 



(2.7) 



and q is the axis ratio. Oblate models have < g < 1, and 
prolate models have g > 1. We write the overall potential 
as a sum $ = + ^ in which the first term is the po- 
tential induced by the stellar density (2.7) and the second 
represents contributions by external components such as a 
central black hole or a dark halo. Different mass densities 
p{m^) require different expressions for the associated poten- 
tial The classical theory of the gravitational potential 
of ellipsoidal bodies (e.g., Chandrasekhar 1969) gives two 
alternative expressions: 



OO oo 



u 
u 



(2.8) 



where G is the gravitational constant, A(m) and U are de- 
fined as 



R^ 



U 



1 + u q^ + u 



and ^0 is the central potential 
2-KGq arcsin e 



/ 



p{m ) dm 



(2.9) 



(2.10) 



with e = -^1 — . In prolate models e is imaginary, but 
arcsin e is real, and equals [q^ — 1)"'"''^ arcsinh(g^ — 1)'"''^. 
The two expressions in equation (2.8) are equivalent only 
when the integral (2.10) converges, and the potential is finite 
everywhere. However, we also need to consider two ways in 
which this integral may diverge. It may diverge only at the 
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lower limit = due to a strong central cusp. In this case 
the potential can be taken as the first expression in equation 
(2.8); it is positive infinite at the centre and vanishes at large 
distances. Alternatively, the integral may diverge only at its 
upper limit = oo. Then the potential can be taken as the 
second expression of equation (2.8), in which $o now is just 
an additive constant. In this case the potential has the finite 
value $0 the centre and becomes negative infinite at large 
distances. When the integral in equation (2.10) diverges at 
both its lower and upper limits, neither expression given in 
equation (2.8) is applicable, since both inner integrals now 
diverge. We must then replace the fixed limits of the inner 
integrals by some interior value of and so take a finite 
part of these divergent integrals. The resulting potential is 
positive infinite at the centre and negative infinite at large 
distances. When convenient, as it is in Section 3.2 below, a 
constant can be added to the potential in all cases. 

The double integration (2.8) can be carried out explic- 
itly in some special cases. More often, only the inner in- 
tegration can be done analytically, and a one-dimensional 
outer integral remains. Some examples are given in Sec- 
tion 3. It is always possible to exchange the order of 
the integration in equation (2.8) to reduce it to a one- 
dimensional integral. While this exchange is simple for B? 
and in the physical range (i.e., both non-negative), it 
must be done more carefully for the wider range of val- 
ues of B? and , which includes complex ones, on which 
our contour integral method operates. For complex val- 
ues of B? and z^ we proceed as follows. Assuming that 
U = -R^/(l + m) + z^ lij^ + m) lies in the region in which 
/o(m^) is analytic, we let ttl' = /(l + a;) + z^/(g^ + a;), 
X £ [u, oo) be the path for the inner complex integral of 
the second expression of equation (2.8). This substitution 
converts the inner integral into one with respect to the real 
variable x. Then an exchange of the order of integration 
followed by a simple integration yields a one-dimensional 
integral. The result is 



also be calculated easily. We find 

oo 

J A(M)(g2 +u) 



oo 

^- = W a(.)(?+.)- ^'(^) + ^- 



(2.12) 



// 2\ —2 /// 2\ —4 

P2 = p{m )q , P22 = p {m )q . 
In order to add a central black hole with mass Mbh, we take 

With these formulas in place, we can choose a pair 
[E, Lz) in the physical domain (Figure 1), and then at each 
quadrature point ^ evaluate B'^ = ^L'l/(^^ — E), solve equa- 
tion (2.6) for z according to the procedure given in Section 
2.1, evaluate equation (2.5), and then compute the contri- 
bution to fe{E, Lz) at the point ^. 

2.3 Velocity profiles 

The observable properties of the two-integral axisymmetric 
models include the line-of-sight velocity moments (e.g., the 
mean streaming velocity {■uios} and the velocity dispersion 
iTios, defined as o-f^j = (■ufoa) — {■wios}^), and the entire VP 
shape. The intrinsic velocity dispersions can be calculated 
conveniently by direct integration of the Jeans equations 
(Hunter 1977; Appendix C). Integration along the line of 
sight can then be done using the expressions given by, e.g., 
Evans & de Zeeuw (1994). The higher order moments can 
be found in a similar way (e.g., Magorrian & Binney 1994). 
Here we discuss only the calculation of the observed VP. 

We let [x, y, z) be Cartesian coordinates with the z-axis 
the symmetry axis of the model. We use [x',y',z') as the 
Cartesian coordinates of an observer, where the z'-axis lies 
along the line of sight, and the x' and y' axes are oriented 
along the major and minor axes of the projected surface 
density of the galaxy. We assume that the galaxy is seen at 
an inclination angle i. Then x = —y'cosi + z'sini, y = x' 
and z = y' sin i -\- z' cos i. 

The normalized VP of an axisymmetric f{E, Lz) galaxy 



^*[B^,z'^) = ^< arcsin e / /)(m^) dm^ + 



oo 

/ 



z' 1 . e 
-, + T":; arcsm — , du > , 

(l+„)2 (g2+„)2j ^2^^^ 



OO 

1 P{V)\ 



oo 

2-KGq f r B^ z^ 



X arcsm e — arcsm 



^/lT 



=) du. 



This allows evaluation of $ by a one-dimensional (numeri- 
cal) quadrature even for complex values of B'^ and z^ . 

The partial derivatives needed in equation (2.5) can now 



YP{vzr,x',y')= ^ J J J f{E,Lz)dv^, dvy, dz', (2.13) 

B>*oo 

where v^i and Vyi are two velocity components on the 
plane of the sky, ■u^/ is the line-of-sight velocity ■uios , and 
S(a;',2/') = J pdz' is the projected surface brightness. We 
employ polar coordinates {v±,ip) in the (■Uj,/ , ■Uj^/ )-velocity 
space, with v^i = v± cos ip and Vyi = v± sin ip. Then 

YP{vz, ■,x',y')=^Jdz' J dvlj f{E, Lz)dv, (2.14) 

zi 

where 

E = nx',y',z')-'^{vl,+vl), 

Lz = -Vz'x' sini (2.15) 

+ v± cos ip \/ {—y' cos i -\- z' sin i)^ + x'"^ cos^ i. 

Here $ = ^{x' , y' , z') is the potential expressed in observer's 
coordinates. The term cos ip should contain a phase shift ipo, 
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but because of the periodicity of the cosine function, and the 
integration over the full range of ip, we can neglect ipo, set 
the range of the (/p-integration to (0, tt), and multiply by two. 
For cases with $<x> = 0, the range of ■u^/ is finite and confined 
by the maximum escape velocity along the line of sight, and 
zi and Z2 are two extreme points along the line of sight de- 



termined by ^{x' , y' , z')- 



0. If*o 



-oo, the range 



of v^i extends from — oo to oo and so does the range of z'. 
Hence zi = — oo and Z2 = oo in this case. We shall generally 
evaluate the triple integral (2.14) by numerical quadrature. 
The even and odd parts of the VP are defined as 

VPe.oK'l 2;', y')=^[YPiv,,;x', y')±YPi-v,,;x' , 2/')]- (2.16) 

According to the above analysis, YP[—v^i]x',y') is given 
by equation (2.14) with only a sign change of the first 
term of the angular momentum in equation (2.15). Since 
f (cos (p)d(p = F{~ cos (p)d(p for an arbitrary function 
F, we can switch the sign of the second term of the an- 
gular momentum in equation (2.15), which then becomes 
— Lz- This relation hence allows us to establish the corre- 
spondence 



2(*-*o 



YPe,o{vz';x 



',y') = ^Jdz' J dvlj f.,o{E,Lz)dv,{2.n) 



The integral (2.14) for the VP simplifies when the term 
—Vzix'smi in expression (2.15) for Lz. vanishes. We show in 
Appendix A that this allows a reduction to a straightforward 
integration over the density distribution itself. This is useful 
for the calculation of three special cases: (i) the VP on the 
minor axis {x' = 0) for arbitrary inclination, (ii) the VP of 
a face-on galaxy (i = 0), and (iii) the density distribution 
of stars which have the systemic velocity, ■u^/ = 0. 



3 SPHEROIDAL MODELS 

We now consider a specific family of models with /)(m^). An 
application of these spheroidal models to the galaxy M32 is 
discussed in Section 4. 



3.1 The (a,/3)-models 

A large number of models with p = p{m^) have been used in 
dynamical studies of galaxies (e.g., Perek 1962; de Zeeuw & 
Pfenniger 1988). We consider the family of models defined 
by 



p{m^) = po 



(3.1) 



where the exponents —3 < a < and /3 < are parameters, 
Po is a reference density, and h a reference length. When 
a = 0, the central density is finite and equal to po- The 
density profile has a central cusp when a < 0. When /3 = 0, 
the models are scale-free. At large radii the density falls off 
proportional to t"'^'^^ . 

When observed at inclination i, these spheroidal models 
have a projected surface density S = S(m'^), given by 



S(m'=) : 



00 



p{m^) dm? 



(3.2) 



where m'^ = x'^ + v'^ /q'^ : so that the isophotes are similar 
concentric ellipses with an observed axis ratio q' which is 
given by 



12 2 • , 2 . 2 • 

q = cos I + q sm i. 



(3.3) 



For the models defined in equation (3.1), S falls off propor- 
tional to (il')""'"^'^"'"'" at large projected radii R' . It has a 
central power-law cusp for a < — 1. When a = —1 this 
cusp is logarithmic. The central projected surface density is 
finite when a > — 1. 

Figure 3 shows a diagram of the (a, /3)-parameter space. 
It is divided into regions by the straight lines a=constant 
and a + 2/3=constant. When a > —2 and a + 2/3 < —2 
(horizontally hatched region), the potential assumes either 
expression of equation (2.8) and its central value is given by 



2'KGqpo\? arcsin e 



B(f+1, 



-/3-1), 



(3.4) 



where B is the beta function. To the left of this region 
(a < —2) the potential is given by the first form of equation 
(2.8), and to the right (a + 2/3 > —2) by the second form. 
When (q!,/3) lies in the vertically hatched region bounded 
by the lines a = — 3 and a + 2/3 = —3, the system has finite 
total mass 



M = 2Tvqpob^B{f + 



3 _a 
2' 2 



-P-l) 



(3.5) 



The line a = — ^ is significant, since a two-integral model 
with a central black hole is physical (/ > 0) only when 
a < — ^ (Section 3.3). 

Some special cases of the family (3.1) have been used 
in dynamical studies before. Along the right boundary of 
Figure 3, models with a = and integer or half integer 
values of /3 are of interest since their potentials are either 
elementary or can be given in terms of special functions (de 
Zeeuw & Pfenniger 1988). The (q!,/3) = (0,-2) models are 
the perfect spheroids (Kuzmin 1956; de Zeeuw 1985) which 
admit three integrals of motion. The scale-free spheroids 
lie along the top boundary (/3 = 0). They are attractive 
candidates for detailed investigation for two reasons. One 
is that their internal dynamics is simpler than general mod- 
els and the other is that they provide good approximations 
to the inner region of cusped models and to the outer re- 
gion of a wide range of spheroidal models. We refer to 
the models with (q!,/3) = ( — 2,0) as the singular isother- 
mal spheroids, since they are the flattened counterparts of 
the well-known singular isothermal sphere. The solid line 
that connects (q!,/3) = ( — 2, —1) to (q!,/3) = (0, —2) indicates 
a set of models that are very similar to ones studied recently 
by Dehnen & Gerhard (1994). Their model densities are like 
equation (3.1), but with (1 + m/fe)^'^ as second term rather 
than {1+m^/b^f. 

The function p(^^,B?) generally can not be given ex- 
plicitly for models with p = p(m?\ This is true also for the 
(a, /3)-family, even in cases where the potential is elemen- 
tary. The two-integral DFs can be found by means of the 
method described in Section 2. The calculations simplify for 
the limiting case of the scale-free spheroids (/3 = 0), which 
we discuss below. 



Axisymmetric galaxy models with central black holes 7 



potentials of the other scale-free spheroids to be 

oo 

(a + 2) J A{u)[ ^ ' \' ^ ' 



Figure 3. The (ck, /3) parameter space that governs the properties 
of the mass models defined in equation (3.1). Scale— free models 
have /3 = 0. Models in the vertically hatched area have finite total 
mass. The central potential is finite in the horizontally hatched 
area. Models with a central black hole have a physical fe(^E , Lz) 
when CK < 0.5, i.e., to the left of the dashed vertical line. The dots 
at (cK,/3) = ( — 2,0) and (0, —2) indicate, respectively, the singular 
isothermal spheroids and the perfect spheroids. The filled square 
at (cK,/3) = (-2,-1) indicates models with a Jaffe (1983) like 
profile. The solid line which connects it to the perfect spheroids 
indicates the set of models that is nearly identical to the family 
studied by Dehnen & Gerhard (1994). The asterisk indicates the 
values of a. and /3 appropriate for the galaxy M32, discussed in 
Section 4. 



3.2 Scale— free spheroids 

The density distribution of the scale-free spheroids can be 
written as 



p = po{B? ^z^ jq) 



.12 



-a ( 

Por (s 



• + q-'' cos' ey'^ , (3.6) 



where R = R/b and z = z/b are dimensionless variables, and 
(f, 6) are scaled polar coordinates defined hy R = f sin 6 and 
z = fcosO. This shows that the density is a product of a 
power of the radius times a function of 5. The total mass of 
these spheroids is infinite. The projected surface density is 
S = So(a;'^ + 2/'V3'^)'^■^"'''^ with q' given in equation (3.3) 
and So =pofe-"-B(|,-f -|)g/g'. 

The potential $o °f equation (2.10) diverges for all 
scale-free spheroids. We therefore replace the fixed limit of 
the inner integral in equation (2.8) by some interior value, 
and make use of the flexibility to add a convenient constant. 
We take the gravitational potential of the singular isother- 
mal spheroid (a = —2) to be 



-■KGqpi 



oo 



du 



In I + z'^ 



(3.7) 



It is oo at the centre and — oo at large radii. It can not be 
expressed in terms of elementary functions, but the associ- 
ated forces can (de Zeeuw & Pfenniger 1988). We take the 



where t) = V j}? and V is defined in equation (2.9). The 
additive constant given by the second term in the bracket 
in expression (3.8) is 



= FoV(« + 2), 
where we have written 
Va = 2TrGqpob^ Ja,q, 



(3.9) 



(3.10) 



and Ja.,q is the auxiliary integral defined in equation (B3). 
The choice of the additive constant $c ensures that the a —> 

— 2 limit of equation (3.8) is equation (3.7). The circular 
velocity ■Uc(-R) in the equatorial plane equals VqR'''^'^^^ , so 
that the value of Vo sets the velocity scale. When —3 < a < 

— 2, the potential is oo at the centre and approaches $c at 
large distances. When a > —2, the potential equals $c at 
the centre and diverges towards — oo at large distances. 

Transformation to the coordinates (f , 6) shows that the 
potentials of the scale-free spheroids, like their densities, are 
all of separable form 



|[lnf= + P_.,,(e)], 



a+2 



(3.11) 



{7=°+=exp[(a + 2)P„,,(e)]-l}, a/-2. 



where Pa.,q{d) is a function of 6 only. It is given explicitly 
in Appendix B. The case a = —2 is a specific example of 
the general set of scale-free potentials considered by Toomre 
(1982). We will work with a scaled potential $ defined as 



r2*/Fo=, a = -2, 

U-*/*=>0, a/ -2, 



(3.12) 



where $c is given in equation (3.9). 

The density (3.6) and its potential can be combined to 
express /5($, R'^) in terms of a function p of a single angular- 
dependent vaiiahle We choose to define this variable as 



R^ exp(*). 



a = -2, 
a / -2. 



(3.13) 



It ranges from on the symmetry axis to 1 in the equatorial 
plane, where Pa,q(7r/2) = 0. We use the definition (3.13) to 
eliminate R dependence from p, and to bring it to the form 



/5(*,-R') 



exp(*). 



(3.14) 



a / -2. 



The definition of the function /o(C) which is introduced here 
is implicit, and is given in equation (B4) of Appendix B. We 
describe there how to evaluate it numerically for complex ^. 

In order to use the contour integral solution (2.2) we 
need the value ^ = $env(-B) where the contour crosses the 
real ^-axis. Use of equations (2.1) gives 



*env(£;): 



2 I Vo\ 



(3.15) 
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It is useful to define two more dimensionless variables, 
namely the scaled energy E and angular momentum Lz, 
as 



bVo' 



2EIVi + 1, 



E : 



(3.16) 



Substitution in the parametric equations (2.1) for the locus 
£ of the circular orbits in the (^E, £2)-plane (Figure 1) then 
shows that the scaled angular momentum Lc of a circular 
orbit with scaled energy E is given by 



Ll{E) 



exp(--B), 

^(a + 4)/(a + 2) 



(3.17) 



OL / -2, 



so that the locus £ can be found explicitly in this case. 

The DFs of the scale-free spheroids are also of special 
form. The contour integral solution (2.2) shows that they 
can be written as 



on , rexp(^), a = -2, 

^ ° a / -2, 

where 



(3.18) 



Figure 4. The ratio fe/ fe{0) as a function of 77^ = 
for scale— free spheroids with axis ratio g = 0.7 and a. ■■ 
— 1.5, and —1. 



-- -2.5, -2, 



L,{E) L.{E)' 



(3.19) 



so that — 1 < 7; < 1. Formulas for calculating values of 
/e(7;^) for the different regimes of a are given in Appendix 
B, together with the elementary expressions for /e(0). These 
formulas are simplified by our choice (3.13) of the variable 
(, which is closely related to the variable rf^ in inversion 
methods (cf. Fricke 1952). 

The velocity moments of the scale-free spheroids can 
be calculated by solution of the Jeans equations, and of the 
higher-order moment equations. We show in Appendix C 
that the second moments and (v^) = (v^) are con- 

nected by a simple relation (eq. [C2]), and that they can be 
expressed in terms of elementary functions when a = —2. 

Equation (3.18) demonstrates that the two-integral DF 
of the scale-free spheroids is a product of a power of energy 
and a function that describes the same relative dependence 
on angular momentum at each energy. This simple form is 
caused by the scale-free nature of the models. The structure 
and dynamics at one radius (energy) are related to those at 
any other radius (energy) by a simple scaling. 

The value of ri indicates the nature of the stellar or- 
bits, from t; = (all orbits with zero angular momentum, 
which are confined to a meridional plane) to t; = ±1 (the 
circular orbits of maximum angular momentum in the equa- 
torial plane). It has been christened the circularity hy Ger- 
hard (1991). Figures 4 and 5 show the ratio /e//e(0) as a 
function of rf^ for different values of a and q. For oblate 
spheroids /e//e(0) is an increasing function of t;^, while for 
prolate spheroids it is a decreasing function oirj^ . The range 
[/e(0), /e(l)] indicates to what extent the high-iz orbits are 
needed to maintain the density distribution of the models. 
At fixed g < 1 this range increases as the density profile 
becomes steeper (a decreases; Figure 4). At fixed a this 
range increases when the flattening of the oblate scale-free 
spheroids increases (g decreases; Figure 5), in agreement 
with results of Dehnen & Gerhard (1994). Flatter models 



Figure 5. Same as Figure 4 but for ol = —2, and q = 0.4, 0.6, 
0.8, 1.0 and 1.2. The scale along the vertical axis is logarithmic. 

require more nearly circular orbits for self-consistent sup- 
port. The importance of the high-iz orbits decreases in 
prolate models: /e(l) drops below zero for sufficiently large 
q, so that prolate two-integral models are physical only for 
a limited range of axis ratios (see below). 

Evans (1993, 1994) constructed a family of self-consis- 
tent axisymmetric models with spheroidal potentials rather 
than spheroidal densities. He took $ oc (il|; + m^)~^^^^ 
with = B? + z^/q'^ and Re and qsi constants. The 
Pb = model has * = -^Vq lii{R% + m^). The density of 
all these models is of the form p{^,R^) = /)o(*) + -R^/02(*), 
and leads to DFs of the form f^{E,Lz) = Fo{E) + LIF2{E), 
where po,P2,Fo, and F2 are elementary functions (powers 
and exponentials). In the limit of zero core radius. Re = 0, 
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these powei — law models are scale-free. The density dis- 
tributions are not spheroidal, become increasingly peanut- 
shaped when qe decreases, and are negative along the z- 
axis when q% < + Pe)- The corresponding DFs are 
similar to expressions (3.18), with a = —2 — Pe, but have 
/e(7;=) = /e(0)[l + ^7;^], where A = A{qE,PE)- Evans' scale- 
free power-law models hence are described by straight lines 
in Figure 4. High-iz orbits are relatively more important 
in our spheroidal scale-free models than in Evans' models. 
This illustrates that in the scale-free models the energy de- 
pendence of /e is determined completely by the slope a of 
the density profile. The 7;-dependence, on the other hand, 
is influenced by a, by the flattening q (or gg), and by the 
shape of the surfaces of constant density. 

The density p(^,B?^ of the scale-free spheroids can 
be expanded in powers of times functions of $ alone. 
Equivalently, the function /o(C) can be expanded in powers 
of C,. Unlike Evans' power-law models, this expansion has 
more than two terms, and, because of the direct relation- 
ship between powers of C, in /o(C) and powers of rf^ in /e (cf. 
Appendix B, eqs. [B8] through [Bll]), the corresponding se- 
ries in Tf^ for /e also has more than two terms. Calculation 
of the successive terms becomes rapidly unwieldy, so that 
evaluation of /e by means of the HQ method, as we have 
done, is more practical. Figure 5 shows that, to first order, 
the scale-free spheroids have /e(7?) ~ /e(0) exp(j47;^), with 
A = A{^q, a) a constant. This approximation is quite accu- 
rate for q in the range between 0.75 and 1.2. It suggests that 
scale-free models with an exact exponential dependence on 
•q^ have nearly spheroidal densities. 

Figure 6 shows the region in the (g, a)-plane where 
the even two-integral DF (3.18) of the scale-free spheroids 
is non-negative. All oblate spheroids of this kind have 
fe{E, Lz) > 0. However, at fixed a there is a maximum 
axis ratio g'max(Q!) > 1 beyond which fe{E, Lz) < for pro- 
late models. This is in harmony with earlier studies of spe- 
cific two-integral axisymmetric models (e.g., Dejonghe & 
de Zeeuw 1988; Batsleer & Dejonghe 1993), and is caused 
by the fact that the t; = orbits needed to reproduce the 
density along the long axis of the model overpopulate the 
density in the equatorial plane when the model becomes suf- 
ficiently elongated. The derived DF corrects this overpop- 
ulation by giving the t; = 1 circular orbits negative weight, 
and hence is unphysical. The range of physical scale-free 
fe{E,Lz) prolate models decreases when a decreases, i.e., 
when the density profile steepens. At fixed q the potential 
becomes more nearly spherical when a decreases, and so do 
the orbital densities, so that the danger of overpopulation 
of the equatorial plane increases. Similar results were found 
by Evans (1994) for the scale-free power-law models. His 
Figure 1 shows a (gB,/3B)-diagram (/3b = —2 — a) which 
can be compared to our Figure 6 (but note that qE is the 
axis ratio of the potential, and not of the density). Physical 
prolate power-law models have a maximum allowed axis ra- 
tio. However, oblate power-law models have fe{E, Lz) > 
only when q% > |(1 + /3b). 

3.3 Small radii: spheroidal cusps and black holes 

For a < 0, the density (3.1) of the (a, /3)-models has a 
power-law cusp near the centre, p Ki po{m/b)'^ , approxi- 
mating the density of the scale-free spheroids. Since the 



Figure 6. The horizontally hatched region indicates the area in 
the (g, ck) plane where the scale— free spheroids have physical DFs 
fe{E,Lz). All oblate (q < 1) models have f^{E,Lz) > 0. At 
fixed CK, physical two— integral prolate models do not exist beyond 
a maximum elongation. The diagonally hatched area indicates 
the more limited region where such spheroids have a physical 
fe{E, Lz) in the presence of a central black hole. 

presence of a central black hole affects significantly the be- 
haviour of the potential near the centre, we discuss cases 
with and without it separately. 

When there is no black hole, the potential can be ap- 
proximated by that of a scale-free spheroid, provided that 
the contribution to the potential from the power-law cusp 
dominates contributions by the rest of the system. This oc- 
curs for systems with parameters (q!,/3) such that a < —2, 
or a > —2 and a + 2/3 > —2, i.e., outside the horizontally 
hatched region in Figure 3. To within an additive constant, 
the potential is given by (3.7) for systems with a = —2, and 
by (3.8) otherwise. For these systems the entire analysis of 
the scale-free spheroids is applicable once the potential $ 
and the energy E are both modified by the additive con- 
stant, so the DF for high energy stars near the centre can 
be readily calculated with the results of Section 3.2 and Ap- 
pendix B. 

Stars outside the central region contribute significantly 
to the finite central potential (3.4) when a > —2 and a + 
2/3 < —2. Then, the correct approximation to the potential 
near the centre is, using the second formula in equation (2.8), 
$ = - + $=usp^ where $="=p and are given in 
equations (3.8) and (3.9), respectively. We now observe that 
p = po{m/b)'^ and $ — $o + form a pair of scale-free 
spheroids with power index a. Hence we only need to replace 
$ and E in the relevant formulas of the previous section, by 
$ — $0 + and E — $S + respectively. For example, 
the dimensionless energy E in equations (3.16) becomes 



E 



2(a + 2) ($; - E) 



When a central black hole is present, its potential over- 
whelms the stellar potential at sufficiently small radii. An 
asymptotic expression for fe{E, Lz) can be calculated as the 
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DF needed to maintain the power-law cusp density in the 
central black hole potential. In most cases the result of in- 
cluding the stellar potential, even in its approximate form, 
is that we are no longer able to obtain a simple, explicit, 
function p{^,R^). Th ere is one exception, which is again 
provided by systems with a > — 2 and a + 2/3 < —2. Then 
the central stellar potential is finite and can be added to the 
black hole potential to approximate the potential at small 
radii by 

GMbh 



+ *0, 



(3.21) 



where $o is the finite central potential (3.4) for the special 
cases mentioned above, and zero otherwise. The first term 
always dominates at sufficiently small radii, but retaining 
the second term provides a more accurate approximation 
with no extra work. Upon solving from equation (3.21) 
and substituting the result into the scale-free cusp density, 
we obtain an approximate, yet explicit, relation 

oc/2 



B 



e R 



V B 



(3.22) 



with B = GMBn/b a reference potential and = 1 — g^. 
This density p{^,R^) is a minor generalization of a compo- 
nent introduced by Dejonghe (1986). HQ give a real one- 
dimensional integral formula for its DF, their eq. (B18), from 
which a factor of 1/2 is missing (Dehnen & Gerhard 1994). 
The DF for the density (3.22) is readily obtained by replac- 
ing E hy E — $0 ill that formula to give 



(2£;-2*o*) 



1/2- 



/ 



-[ic-".)»l[s^-.v]-'-«| 

poq /' -fc - Wq ^ ' F / 2 2n 



(3.23) 



£3/2 V B 

where the function is an elementary real integral, which 
can be written in terms of a generalized hypergeometric 
function (Dejonghe 1986) 

1 r(i-a) 

(3.24) 



(27r)3/2 r(-a- 1/2) 



X 3-F2( — ,1 - J.-j;-" - 2' 2^^ ^ )■ 

The DF (3.23) is again of a separable form: the product of a 
power of [E — ^q)/ B and the function whose argument is 
e=7,= = eni/LliE), where LUE) = (GMbh)V[2(-B - *o*)] 
is the square of the angular momentum of a circular orbit 
of energy E. The function depends only on the power 
a of the density cusp, and not on the axis ratio q, though 
this ratio is part of the argument of /„. For integer a the 
generalized hypergeometric function in the definition of 
reduces to an ordinary hypergeometric function, which can 
be written in terms of elementary functions. Table 1 lists 
fa for the cases a = — 1, — 2 and —3, which have very simple 
forms. The case a = —4 is given in Appendix B of Dehnen & 
Gerhard (1994). For integer a < —4 the expressions rapidly 
become lengthy. 

The value of /ci(0) is the elementary factor in front of 
3F2 in equation (3.24), so that /e(£^, 0) is elementary. For 



Table 1. Some special cases of fa. 
« /c(e=7,= ) 

2^/2^^^ (1-^2^2)2 



(2-t- 77^ L— e2 7^2_j.3g^ arcsin 677 



(l_e2^2)5/2 



2\/2 l + 3e^77^ 
71-2 (l_e2,,2)3 



spherical cusps [q = 1, e = 0), fe{E, 0) gives the isotropic 
flE) for large energy E. It agrees with the expression given 
by Tremaine et al. (1994) for Dehnen's (1993) family of 
cusped spherical models (but they neglect the term $o)- 
For oblate cusps is a monotonically increasing function 
of its argument, and the DF (3.23) is least for stars with 
zero angular momentum (t; = 0). The extreme value /a(e^) 
becomes large as e grows and the isodensity contours flatten. 
For /e(£^, 0) to be non-negative, we must have a < —1/2. 
This is the same constraint on the slope of the cusp density 
profile as in the spherical case (Tremaine et al. 1994). In 
prolate cusps = 1 — < 0, and the DF (3.23) is least for 
stars on circular orbits (t; = 1). It is, of course, unphysical 
for a > —1/2, but for smaller a it is physical only when 
q is less than a certain maximum value. The region in the 
[q, a)-plane where spheroidal cusps have physical fe{E, Lz) 
in the presence of a central black hole is shown in Figure 
6. For prolate models this allowed region is more limited 
than the area occupied by physical self-consistent f{E, Lz) 
scale-free spheroids. Near the black hole the stars experi- 
ence a spherical potential. As we have seen in Section 3.2, 
this increases the danger of overpopulating the density in 
the equatorial plane of a prolate model, so that at fixed a 
the physical range of g > 1 is smaller. The potential of 
the self-consistent spheroids becomes more nearly spherical 
(and Keplerian) when a decreases towards —3 at fixed q, so 
that the shrinking of the allowed range of q caused by the 
inclusion of the black hole is less. 

3.4 Large radii: power— law halos 

At large radii the density approaches p = po{m/b)°''^'^^ , 
which is again scale-free. We need to distinguish systems 
of finite mass (a + 2/3 < —3; the vertically hatched area in 
Figure 3) and those of infinite mass (a + 2/3 > —3). We 
ignore the intermediate a + 2/3 = —3 case in which the total 
mass is marginally infinite. For the former, both the stellar 
potential and that of any central black hole are Keplerian 
and both should be included in an approximate potential, 
unless Mbh <C M. The DF is again given by formula (3.23) 
with $0 = 0, a replaced by a + 2/3, and the total stellar 
mass M (3.5) added to Mbh in the definition of B. 

In a self-consistent system of infinite mass, the stellar 
potential at large radii dominates any Keplerian potential, 
and so is insensitive to any central black hole. Therefore the 
potential approximation for cases with or without a central 
black hole is the same. Hence, all systems with a + 2/3 > —3 
are approximated at large radii by scale-free spheroids of the 
type discussed in section 3.2; the approximate DF fe{E, Lz) 
for stars at large radii [E —> $<x>) can be obtained once we 
replace a by a + 2/3 and modify both $ and E by additive 
constants if necessary. 
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flattened 



We have also investigated the case of a 
spheroidal density p = po^m/b)" , with 
in the non-self-consistent scale-free power-law potential 
defined as (cf. Evans 1994) 

( Inmd, 



0, 



(3.25) 



where md 



i^-l], 7/0, 
TTid/c, with = B? + so that qd is 



the axis ratio of the spheroidal equipotentials and c is a 
reference length. When qd = 1 and 7 = — 1 this reduces 
to the case discussed in Section 3.3 of a spheroidal density 
po{m/b')'^ in the Kepler potential. Proceeding as in that 
case, we again find an explicit expression for p{^d, R^) which 
is another generalization of the Dejonghe (1986) component 
(cf. eq. [3.22]). When 7 / 0, it is given by 



1/7 



2 52 

SdRd 



7*d 



2/71 



(3.26) 



where = 1 — q'^/q'd and Rd = R/c. When 7 = 0, so that 
the potential is logarithmic, we find 

p{<id,R') = Po(y^] exp(-a*d/Fo')x 



cqd 
V bq 

[1 -e^^^exp(2*d/Fo= 



1 "/2 
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Application of the HQ method shows that the corresponding 
/e is given by 



'exp(-^), 7 = 0, 



7/0, 



where E is defined as 



2M. 4. 1 



E 



2 

7 + 2 



(1-^), 7/0. 



(3.28) 



(3.29) 



The /a, 7(6^7?^) can be evaluated in the same way as the 
function /e of equation (3.18), as detailed in Appendix B. 
These DFs are useful for understanding the behaviour of 
two-integral models in the limit of large radii, where the 
density falls off as a power-law, and the potential may be 
dominated by a dark halo. We remark that changing the 
axis ratios q and qd as well as the normalisation lengths b 
and c while keeping both q/qd and b/c constant leaves the 
DF invariant. 



4 APPLICATION TO M32 

To illustrate our technique we have used it to model the 
galaxy M32, which is believed to harbour a massive black 
hole in its nucleus (Tonry 1987; Dressier & Richstone 1988; 
Richstone, Bower & Dressier 1990). HST observations have 
revealed the presence of a central surface brightness cusp 
(Lauer et al. 1992). High spatial resolution (ground-based) 
kinematical and VP measurements along several position 
angles are available from van der Marel et al. (1994a, here- 
after vdM94a). Axisymmetric f{E, Lz) models were used by 
van der Marel et al. (1994b, hereafter vdM94b) to interpret 



Figure 7. Contour plot of the even part fe{E,Lz) of the DF 
of our model for M32, which has a central black hole of mass 
1.8 X 10^ M0, as discussed in the text. The quantities along the 
axes are 77= = Ll/Ll{E), and the energy E in units of , which 
is the central potential in the absence of a black hole. The allowed 
energy range is £/ G [0,oo), but only the range [0.1, 2] is shown. 
At lower and higher energies the DF has reached its asymptotic 
behaviour as dictated by the scale— free approximations. Adjacent 
solid contours are a factor 0.17 apart, the highest 'contour' being 
at the maximum of the DF, which is indicated by the solid dot. 
The ten dashed contours are a factor (0.17)°'^ apart, the highest 
contour again being at the maximum of the DF. 

these observations. The modelling consisted of: (i) use of 
Evans' (1994) power-law model DFs without a central black 
hole; and (ii) calculation of the first three moments of the 
VP for an (a, /3)-model with a black hole, by solution of the 
moment equations of the coUisionless Boltzmann equation. 
A remarkably good fit was obtained with a black hole of 
mass Mbh ~ 1.8 x 10^ M©, but the actual DF of the model 
could not be calculated. With the technique presented here 
we can calculate the entire DF (Section 4.1). This in turn 
yields the full VP shapes, and hence allows a more accurate 
comparison to the data of vdM94a (Section 4.2). It also al- 
lows a detailed study of the expected VP shapes for the high 
spatial resolution spectroscopic observations that are soon 
to be expected with the HST (Section 4.3). 



4.1 The f{E, Lz) distribution function for M32 

The observed surface brightness distribution of M32 can 
be well fitted with an (a, /3)-model with parameters a = 
-1.435, P = -0.423, b = 0.55", po = ioTyMo/Lo ^,, 
jo = 0.470 X 10^{q' /q)LQ y pc"^ and q' = 0.73. This leaves 
three free parameters that can be chosen to optimize the 
fit to the kinematical observations: the inclination i (which 
enters in the relation [3.3] between q and q'), the average 
T^-band mass-to-light ratio of the stellar population in 
solar units (assumed to be independent of radius), and the 
black hole mass Mbh. We restrict ourselves here to the 
model that is edge-on (i = 90°) and has = 2.51 and 
Mbh = 1.8 x 10^ Mq, based on the results of vdM94b. We 
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do not discuss the detailed dependence of the model pre- 
dictions on the parameters i, and Mbh, since we do 
not expect the results of our models to change the discus- 
sion given previously by vdM94b, on the basis of their more 
approximate modelling. 

Figure 7 shows a contour diagram of the even part 
fe{E, Lz) of the DF for our model for M32, obtained with 
the technique described in Sections 2 and 3. The DF is pos- 
itive for all physical values of [E, Lz)- The contours slope 
gently upwards in the right half in accordance with the weak 
dependence of /e on energy predicted by the scale-free ap- 
proximation (3.23) for our model with a = —1.435 (which 
is close to —1.5, the value for which the energy dependence 
in eq. [3.23] vanishes). The contours slope sharply down- 
wards in the left half in accordance with the strong depen- 
dence of /e on energy predicted by the scale-free approxi- 
mation (3.18) for the appropriate a value of —2.281. Fig- 
ure 8, which shows the dependence of /e on the energy E, for 
both Lz = and Lz = Lci^E), confirms the accuracy of the 
scale-free approximations for the limit of low and high en- 
ergies. Figure 9 shows a contour plot of fe{E, Lz)/ feiEjO). 
The nearly horizontal contours indicate that the /e for our 
M32 model is close to being a separable function of E and 
Tj'^ = L^/ L^[E) at all energies, with relatively small dis- 
crepancies in the transition region 0.5 < E < 1.1 between 
low and high energies. The similarity of the behaviour of 
fs{E,Lz)/fs{E,0) as a function of t;^ in the low and the high 
energy limit is further illustrated in Figure 10, which shows 
the asymptotic behaviour obtained from the scale-free ap- 
proximations. For comparison, the predicted behaviour at 
high energies for the same model without the central black 
hole is also shown. The dependence of fe{E,Lz)/fe{E,0) 
as function of t;^ is much steeper at high energies when the 
black hole is present, because then the stars close to the 
centre orbit in a spherical rather than a flattened poten- 
tial. This causes the density contributed by stars of the 
same [E, Lz) to be more nearly round than in the case of a 
flattened potential (Section 3.4). In order to reproduce the 
same flattened density distribution, the number of stars on 
high-iz orbits must therefore increase relative to the case 
where the potential is flattened. 

Apparently, the mass density slope in the inner parts of 
M32 and the presence of the black hole 'conspire' to produce 
nearly the same dependence of fe{E, Lz)/fe{E, 0) on rf^ at 
high energies, as at low energies, where the dependence is 
governed by the mass density slope in the outer parts. It is 
not clear whether this is a mere coincidence, or is the result 
of stellar dynamical processes which have operated in M32, 
possibly caused by the presence of the central black hole. 
Such processes would then have to be capable of removing 
any dependence of the DF on a third integral of motion 
and would have to drive the DF to a product form, all in 
less than the Hubble time. More quantitative theoretical 
work, e.g., by study of the adiabatic growth of central black 
holes in stellar systems (Young 1980; Quinlan, Hernquist & 
Sigurdsson 1994), is clearly needed in this area. 



4.2 Comparison to ground— based kinematical data 

To compare the kinematical predictions of our model to the 
available data we must specify also the odd part fo of the 



Figure 8. Behaviour of fe{E,Lz) for our model for M32 dis- 
cussed in the text, as function of energy E, for 77^ = [Lz = 0; 
solid line) and for 77^ = 1 [Lz = Lc{E); dashed line). The dot- 
ted lines show the asymptotic slopes at low and high energies 
expected from the scale— free approximations. These slopes agree 
well with the model slopes. The actual values predicted by the 
scale— free approximations also agree well with the model predic- 
tions (for all 77^ ) , but are not shown here for the purpose of clarity. 



Figure 9. Contour plot of fs{E,Lz)/ fs{E,0) for our model for 
M32 discussed in the text. The quantities along the axes are as 
in Figure 7. Adjacent contours in the figure differ by a factor 0.8. 
The highest contours lie at the top of the plot. 



DF. We choose the parametrization 
f.{E,Lz) = {2F 



tanh[a7?/2] 
1) .„„v.r„/m f-i^' 



tanh[a/2] 



Lz), 



(4.1) 



where as before 7; = Lz/ Lc{E), and < f < 1 and a > 
are free parameters. This is a modified version of a func- 
tional form derived by Dejonghe (1986) based on maximum 
entropy arguments. With this choice for fo, the total DF 
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Figure 10. The behaviour of f^{E,L;,)/ f^{E,0) for our model 
for M32 discussed in the text, as a function of 77^ = L^/L^[E), in 
the Umit of low energies (solid curve) and high energies (dashed 
curve). The dotted curve shows what the behaviour in the limit 
of high energies would have been without the central black hole. 

/ = /e + /o is positive whenever /e is. We adopt f = 1 
and a = 5.5, based on the results of vdM94b for Evans' 
power-law models. The DF of our model for M32 is now 
specified completely, and the model VPs can be calculated 
as described in Section 2.3. Figure 11 compares the model 
predictions to the data along five different slit positions pre- 
sented by vdM94a. Each VP is characterized by six parame- 
ters: the mean V and dispersion a of the best-fitting Gaus- 
sian to the VP, and the Gauss-Hermite moments hs, . . . , he, 
defined as in van der Marel & Franx (1993). Our model 
predictions take the spatial binning and seeing point spread 
function (PSF) convolution of the observations into account, 
as described in Appendix D. Small (< O.l") offsets of the 
slit from the galaxy centre due to differential atmospheric 
refraction were also modelled. The results in Figure 11 con- 
firm the main conclusions of the modelling by vdM94b. The 
model fits the data remarkably well, much better than one 
would have expected a priori. The observed steep central 
gradient in the mean velocity and the observed central peak 
in the velocity dispersion are both reproduced (owing to 
the presence of the central black hole in the model). The 
observed Gauss-Hermite coefficients are fitted up to a rms 
residual of only ~ 0.02, indicating that the dynamical struc- 
ture of M32 is most likely very close to that of an f{E, Lz) 
model. Nonetheless, some minor discrepancies between the 
observations and the model predictions remain, most likely 
indicating a (slight) dependence of the DF on a third inte- 
gral. 

First, the small discrepancies between the observed and 
the predicted Gauss-Hermite coefficients outside the central 
arcsec on the major axis, are in the sense that the even part 
of the observed major axis VPs is slightly more flat-topped 
than predicted. This might indicate that M32 has a veloc- 
ity distribution with (v^) > (vg) > (I'r)- The discrepancy 
is hardly significant, however, and is in fact smaller than 
that inferred by vdM94b. This we have found to be due to 



the fact that the (a, /3)-model used here has an asymptotic 
mass-density slope that is slightly steeper than that of the 
power-law model employed by vdM94b (p oc m~^'^^ versus 
p oc m~^'^). It is not a consequence of the fact that the iso- 
density contours of the (a,/3)- and power-law models have 
slightly different shapes. 

Secondly, Figure 11 shows that the predicted amplitude 
of the mean line-of-sight motion V along the intermediate 
axes is slightly too high. This is consistent with the find- 
ings of vdM94b, who argued that to obtain a good fit on 
both axes, one must invoke an odd part of the DF that 
depends also on a third integral. To test this conclusion, 
we attempted to solve the inverse problem. We adopted a 
streaming velocity field defined by (v^)^ = k(^(v^) — (v^)) 
with k = 1.25, which provides a good fit to {■uios} on both the 
major and the intermediate axes (vdM94b). We then used 
the HQ method to obtain the unique fo{E,Lz) consistent 
with this streaming velocity field. We found that the result- 
ing model is unphysical, because the total DF / = /e + /o is 
not positive for all physically accessible [E, Lz)- It thus ap- 
pears indeed, that an odd part of the form fo = fo{E, Lz, I3) 
is required to fit the amplitude of the mean streaming mo- 
tions along all slit positions simultaneously. 

Thirdly, even with the inclusion of a central black hole 
there remains a discrepancy between the observed and the 
predicted velocity dispersions near the centre. Especially 
on the minor axis, the observed central peak in the velocity 
dispersion is steeper than that predicted by the model. This 
hints at a dependence of the even part of the DF on a third 
integral. 

The most important conclusion from the f{E, Lz) mod- 
elling of M32, however, is that, aside from the minor discrep- 
ancies discussed above, the accuracy with which an f{E, Lz) 
model can fit the data is quite remarkable. A similar conclu- 
sion was reached independently by Dehnen (1995), who used 
a Richardson-Lucy algorithm (Newton & Binney 1987) to 
construct f{E, Lz) models for M32. Our results are not very 
sensitive to the assumed inclination angle (vdM94), which is 
not well constrained by the data. It is also not very sensitive 
to the precise value of the slope of the density profile inside 
0.3", which might differ slightly from the value adopted here 
(Lauer et al. 1992). With a slightly steeper slope an equally 
good fit is obtained with a slightly less massive black hole, 
and vice-versa. 

4.3 Predictions for observations with HST 

The fact that an f{E, Lz) model with a central black hole 
can provide such a good fit to ground-based kinematical 
and VP data does not necessarily imply that M32 must 
have a central black hole. To date it has not been convinc- 
ingly demonstrated that three-integral axisymmetric mod- 
els without a central black hole cannot also fit the same 
ground-based data. High spatial resolution data from the 
refurbished HST should provide more definite evidence ei- 
ther for or against the presence of a central black hole. In 
this section we discuss the kinematical and VP predictions 
of our model for M32, for observations through the small 
apertures available on the HST. This yields definite predic- 
tions for the signatures of a central black hole that one might 
expect to observe with the HST. 

We discuss the normalization 7, the mean V and the 
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Figure 11. The data points are the observed kinematics and VP parameters for M32 as a function of the projected radius R' along five 
different sUt position angles, as presented by vdM94a. From top to bottom: the mean and dispersion of the best— fitting Gaussian to the 
VP, and the Gauss— Hermite coefficients /13 , . . . , /ig . The curves show the predictions of our model for M32, which has a 1.8 X 10^ 
central black hole, taking into account the seeing convolution and spatial binning of the observations. 
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Figure 12. Predicted Icinematics and VP parameters for our 
model for M32, wliicli lias a 1.8 X 10^ M© central black hole, 
for observations through a circular aperture placed on the galaxy 
centre. Solid curves show as functions of the aperture diameter Z): 
the normalization 7 and dispersion a of the best— fitting Gaussian 
to the VP, and the Gauss— Hermite moments /14 and he . The 
short— dashed curves in the left panels show the true normalization 
and dispersion of the VP. 

dispersion a of the best-fitting Gaussian to the VP, as well 
as the Gauss-Hermite moments up to order 6. In a real ob- 
servational situation a galaxy spectrum is modelled as the 
convolution of a (stellar) template spectrum and a broaden- 
ing function. The ratio of the equivalent width of the absorp- 
tion lines in the galaxy spectrum to those in the template 
spectrum is called the 'line strength'. This line strength is 
unknown, and has to be estimated from the data. If this 
true line strength is one, then our parameter 7 is the esti- 
mate of the 'line strength' that one would expect to obtain 
by fitting a Gaussian broadening function to the data. 

4-3.1 Predictions for centred apertures 

We first calculated the predicted VPs for observations 
through a circular aperture placed on the galaxy centre. Fig- 
ure 12 shows the predicted (7, it, hi, he) as functions of the 
aperture diameter D [the quantities (V, hs, he,) are zero since 
the central VPs are symmetric]. The true line strength (as- 
sumed to be = 1) and dispersion of the VP are shown also. 

The true velocity dispersion of the VP increases with 
decreasing D, roughly as ir^ Ri ci -\- C2/D, where ci and C2 
are constants. For D > 0.5", the predicted VP is close to 
Gaussian. For smaller diameters the VP wings are more ex- 
tended than those of a Gaussian (see also Figure 15 below). 
This is due to the stars that orbit close to the hole at very 
high velocities, and is quantified by the increasingly non-zero 
values of hi and he . The non-Gaussian wings of the VP con- 
tribute significantly to the normalization and dispersion of 



Figure 13. Predicted kinematics and VP parameters for our 
model for M32, which has a 1.8 X 10^ M© central black hole, 
for observations through a 0.09" X 0.09" square aperture, placed 
along the major axis at a distance from the galaxy centre. 
Solid curves show as functions of R': the normalization 7, mean 
V and dispersion a of the best— fitting Gaussian to the VP, and the 
Gauss-Hermite moments /13 , . . . , /ig ■ The short-dashed curves in 
the left panels show the true normalization, mean and dispersion 
of the VP. The ground-based major axis V and a observed by 
vdM94a are also shown for comparison (dots). The long-dashed 
curve shows the model fit to these data. The results illustrate the 
major improvement to be expected with the HST. 

the VP. A Gaussian fit is insensitive to the wings of the VP, 
and hence underestimates both the true line strength and 
the true velocity dispersion. 

Our predictions for M32 are qualitatively similar to 
those of van der Marel (1994d), who discussed the expected 
kinematics and VP shapes for centred aperture observations 
of the galaxy M87, using a spherical model with a 5 x 10^ Mq 
central black hole. Nonetheless, there are a few noticeable 
differences. The stars that are influenced significantly by a 
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Figure 14. As Figure 13, but now for observations through a 
circular aperture of diameter D = 0.26", placed along the major 
axis at a distance from the galaxy centre. 

central black hole in a stellar system reside in a sphere with 
radius of order t-bh = 2GMBH/37rCT^, where F is a 'typical' 
velocity dispersion outside the region influenced by the black 
hole. With this definition, the projected velocity dispersion 
IT of a singular isothermal sphere with a massive black hole 
satisfies: = + ('•bh/-R')] (Tremaine et al. 1994). In 
arcseconds on the sky: 

where d is the distance to the galaxy. This radius is approxi- 
mately four times smaller for M32 (Mbh w I.SxIO^Mq, F w 
55kms-\ d w 0.7 Mpc) than for M87 (Mbh w 5 x IO^Mq, 
W w 300 km s"'', d Ri 16 Mpc). So to detect similar devia- 
tions from a Gaussian for both galaxies, M32 must be ob- 
served with a four times smaller aperture than M87. Con- 
versely, for a fixed aperture size, the expected VP for M32 



Figure 15. The two solid curves are the VPs predicted by our 
model for M32, which has a 1.8 X 10^ M© central black hole, 
for observations through a 0.09" X 0.09" square aperture: (i) 
placed on the galaxy centre (i?' = O"); and (ii) placed along the 
major axis at R' = 0.1". Both VPs are normalized. The two 
dotted curves are the best— fitting Gaussians to these VPs. The 
arrows indicate the central escape velocity ±-^2^^, due to the 
gravitational potential of the stars. If there were no black hole in 
the model, no stars would be observed beyond this velocity. 

is much closer to a Gaussian than that for M87. 

From an observational point of view, there are two main 
differences between M32 and M87. First, M32 has a much 
smaller velocity dispersion. So while for M87 the contin- 
uum subtraction in the spectral analysis is a serious prob- 
lem (van der Marel 1994c, d), this is not expected to be 
the case for M32. On the other hand, for M32 the lim- 
ited instrumental resolution will be a complicating factor. 
The Faint Object Spectrograph (FOS) aboard the HST has 
finstr ~ 100 km s"'', of the same order as the stellar velocity 
dispersion. The second important difference between M32 
and M87 is that M32 has a much higher surface brightness. 
For observations of the M87 centre with a. D = 0.26" aper- 
ture, exposure times of > 10 hours are required to obtain a 
sufficient signal-to-noise ratio for a useful VP analysis. For 
M32, not more than ~ 15 minutes are required. 

4-3.2 Predictions for apertures placed along the major axis 

To obtain constraints on the rotational properties of M32 it 
will be useful to obtain HST aperture observations at various 
distances along the major axis. We therefore calculated the 
predicted VPs of our model for M32, for observations with: 
(i) a 0.09" X 0.09" square aperture; and (ii) a D = 0.26" 
circular aperture. The first is the size of the smallest aper- 
ture available on the HST/FOS, the second is the size of 
the smallest circular aperture available on the HST/FOS. 
Figures 13 and 14 show the predicted (7, V, it, /13, . . . , fee) as 
functions of the galactocentric distance R' on the sky. As 
in Figure 12, the true normalization, mean and dispersion 
of the VP are shown for comparison, as well as the observed 
and predicted ground-based V and it (from Figure 11). Fig- 
ure 15 shows the predicted VPs with the square aperture for 
R' = 0" and R' = O.l", together with the best Gaussian fits 
to these VPs. 

The velocity dispersion one would expect to measure 
in the centre [R' = O") by fitting a Gaussian VP to the 
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data, is 127 km s"'' for the square aperture and 105 km s"'' 
for the circular aperture. This is significantly larger than 
the central velocity dispersion of ~ 86kms~'' obtained from 
ground-based data (see Figure 11). Figure 15 shows that the 
broad wings of the central VP provide a strong signature of 
the black hole. The arrows in the figure indicate the central 
escape velocity ^2$o due to the gravitational potential of 
the stars, which for our model is 298 km s"'^. In the absence 
of the central black hole no stars would be observed beyond 
this velocity. 

Outside the centre [R' ^ 0) the predicted VPs are 
asymmetric, with a tail away from the direction of rotation 
(y and hs of opposite sign). This is evident in Figure 15. 
Similar VPs are observed from the ground (see Figure 11). 
The mean of the best-fitting Gaussian overestimates the 
true mean velocity by about 15%, as a result of the VP asym- 
metry. The mean streaming curves in Figures 13 and 14 have 
similar shapes. They rise almost linearly out to a character- 
istic radius which is of the same order as the aperture size, 
and then remain flat at ~ 50 km . The circular velocity 
of the model has a pronounced Keplerian [vc oc r~^^^') in- 
crease close to the black hole. However, only a very minor 
increase is seen in the predicted mean streaming curve in 
Figure 13, and no increase at all is seen in Figure 14. The 
reason for this is that even the smallest aperture available 
with the HST/FOS is not much smaller than the radius t-bh 
defined in equation (4.2), which for M32 is ~ 0.1", 

For R' = O.l", the mean velocity of the best-fitting 
Gaussian to the VP is 50 km for the square aperture 
and 35 km for the circular aperture. If, this close to the 
centre, mean streaming velocities of this order were actually 
measured with HST, it would probably provide a strong ar- 
gument against models without a central black hole. Such 
models require a large amount of radial motion close to 
the hole to account for the high central velocity dispersion 
(vdM94b), and in such models the maximum possible mean 
streaming is limited (Richstone et al. 1990). 

5 SUMMARY AND CONCLUSIONS 

The contour integral method of Hunter & Qian (1993) can be 
used to calculate the even part f,{E, Ll) of the DF f{E, L,) 
for smooth axisymmetric densities p(^R?,z^) in a potential 
$(il^,z^). Unlike previous methods, the HQ method is ap- 
plicable in cases where p as a function of $ and R? — de- 
noted here by , R?) — is not known explicitly, and this 
key property finally allows the construction of large classes of 
realistic axisymmetric galaxy models. We have shown how 
this can be accomplished for the family of classical spheroids, 
in which the density distribution is stratified on similar con- 
centric oblate or prolate spheroids with constant axis ratio 
and has an arbitrary radial profile. In projection, these mod- 
els have concentric elliptic isophotes with constant elliptic- 
ity, and no isophote twist. The "density" p(^,R?^ of these 
models is generally only known implicitly. The HQ method 
requires evaluation of p at complex values of R? and and 
we have described in Section 2 how this analytic continua- 
tion can be done numerically. It is then straightforward to 
evaluate the contour integral for fe{E, Lz)- 

Our procedure for the calculation of f{E, Lz) applies 
not only to a single spheroidal component, but also to any 
combination of them with different axis ratios and density 



profiles. In particular, it can be applied to the sums of Gaus- 
sian density distributions that have been used to represent 
rather complicated axisymmetric models of realistic galaxies 
(Monnet, Bacon & Emsellem 1992; Emsellem et al. 1994). 
In a future paper we shall use our method to provide two- 
integral DFs for such models and to obtain further insight 
into their structure and dynamics. It remains to be seen 
whether f{E, Lz) can be calculated in an analogous way for 
arbitrary smooth axisymmetric densities that are not of the 
form p = /)(m^), for example by direct use of the Poisson 
integral for the potential (Binney & Tremaine 1987, eq. [2- 
3]). 

In Section 3 we considered a specific set of classical 
spheroids. These (a,/3)-models have arbitrary axis ratio, 
and the slopes of the density profile in the inner and outer 
regions can be chosen independently. The (a, /3)-family con- 
tains many popular axisymmetric models as special cases, 
including the scale-free spheroids in which the density pro- 
file is a pure power law. The "density" /5($,il^) of these 
models is not known explicitly, but it has a simple form, 
and application of the HQ method is straightforward. We 
have calculated the resulting DFs for a variety of axis ratios 
and density profile slopes. 

The scale-free spheroids can be compared to the scale- 
free power-law models of Evans (1993, 1994), for which the 
potential rather than the density is stratified on similar con- 
centric spheroids. Evans' models have elementary fe{E, Lz), 
which lead to elementary and explicit expressions for the ob- 
servables (Evans & de Zeeuw 1994). However, their density 
distributions can deviate strongly from a spheroidal shape, 
and may be peanut-shaped (even though this is less evident 
in projection). The non-spheroidal shape of these models 
is reflected in a linear dependence of ME,Lz) on Ll. By 
contrast, the DF of the scale-free spheroids has the same 
energy dependence as the Evans models — which is fixed 
by the slope of the pure power-law density profile — but 
the dependence on angular momentum is stronger, so that 
the high-iz orbits are more heavily populated. The ad- 
vantage of the scale-free spheroids presented here is that 
they have exactly spheroidal density distributions, but this 
pleasing property comes at a price: the DF and the observ- 
ables are not elementary functions, and require numerical 
integrations, which are however straightforward. The scale- 
free spheroids can be used to approximate the behaviour of 
the general (a, /3)-models at small and large radii. Their 
simpler structure speeds up the calculation of f{E, Lz), and 
hence allows an efficient investigation of parameter space. 
We have determined in Sections 3.3 and 3.4 the flattenings 
and density profile slopes for which oblate and prolate cusps 
have physical, i.e., non-negative, two-integral DFs. We have 
shown that the physical set of two-integral prolate models 
is limited in axis ratio, and shrinks even further when a 
black hole is included in the potential. We also extended 
the computation of fe{E,Lz) of the self-consistent scale- 
free spheroids to the case of scale-free spheroidal densities 
embedded in Evans' power-law potentials of arbitrary flat- 
tening and radial profile. These DFs allow a systematic in- 
vestigation of the effect of a dark halo on the observed VPs 
in flattened elliptical galaxies, and hence should be useful 
for the analysis of kinematic measurements at large radii. 

High-resolution ground-based kinematic measurements 
for the galaxy M32 by vdM94a were interpreted by vdM94b 
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in terms of an (a,/3)-model with f(^E,Lz) and a 1.8x10^ Mq 
central black hole, for which they solved the second and third 
order moment equations of the coUisionless Boltzmann equa- 
tion. In Section 4 we have computed the exact two-integral 
DF for this model. We used the HQ method to calculate /e, 
chose a simple functional form for fo that fits the observed 
mean streaming velocities {■uios}, and computed the expected 
VPs for edge-on observation, taking into account the seeing 
convolution and spatial binning of the observations. The 
results confirm that this f{E, Lz) model provides a truly re- 
markable fit to the available data. In addition, it turns out 
to have a remarkable property: it is close to being a sim- 
ple product of a function of energy times a function of the 
circularity parameter ri = Lz/Lc, with Lc the angular mo- 
mentum of the circular orbit with energy E. Quantitative 
theoretical work is needed to determine whether this result 
has any important physical significance. 

The success of a two-integral model with a central black 
hole is no proof that M32 indeed contains such a black hole, 
as we have not demonstrated that a three-integral axisym- 
metric (or triaxial, see Emsellem et al. 1993) model without 
a black hole can be ruled out. Spectroscopic observations 
with the high spatial resolution of the HST should provide 
more definite evidence either for or against the presence of a 
central black hole. We have used our model for M32 to pre- 
dict what HST should reveal. We calculated the expected 
VPs for spectroscopic observations with the smallest rect- 
angular (0.09" X 0.09") and circular [D = 0.26") apertures 
available on the HST/FOS. The predicted central Gaussian 
velocity dispersion is 127 km with the former, and 105 
km with the latter aperture. It is not expected that 
one will be able to measure the expected Keplerian rise of 
{■uios} close to the black hole, because its radius of influence 
is only ~ O.l". When measured with the available small 
apertures, the predicted mean streaming motions along the 
major axis are nearly constant at ~ 50 km s"'", down to 
O.l". If such mean streaming motions are indeed measured 
at O.l" from the centre of M32, then it will be very hard to 
argue for models without a central black hole. Such models 
require a strongly radially anisotropic velocity distribution 
near the centre in order to account for the observed large 
central velocity dispersion, and hence cannot support large 
mean streaming motions. 

It has finally become practical to calculate f{E, Lz) for 
realistic axisymmetric galaxy models. We have shown here 
that one way to do this is to use the HQ method. Other pos- 
sibilities include the series expansion method of Dehnen & 
Gerhard (1994) and the grid-based quadratic programming 
technique of Kuijken (1995). As a result, two-integral ax- 
isymmetric models can now replace spherical isotropic mod- 
els as the standard theoretical templates for a zeroth order 
comparison to the high quality kinematic observations of 
flattened elliptical galaxies that are available. The case of 
M32 shows that f{E, Lz) modelling may already provide a 
remarkable fit to some galaxies, but it is well-known from 
modelling based on the Jeans equations that this must be 
the exception rather than the rule. Application of these im- 
proved modelling techniques to elliptical galaxies with more 
internal structure, such as those with embedded discs, will 
be quite rewarding. 
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APPENDIX A: THREE SPECIAL VELOCITY 
PROFILES 

When the term —v^ix' sini in equation (2.15) vanishes, the 
innermost integral of expression (2.14) for the VP spans an 
interval of Lz which is symmetric about Lz = 0. Conse- 
quently only /e contributes to it. Changing the integration 
variable to i^, the two inner integrals can be written as 
follows 



« ) 



I"" I 



U{E,Lz)dLl 



/i| V2s=(^ -E)- LI 



(Al) 



where 



1 2 



= {—y' cos i + z' sin i)'' + x cos" i. 



(A2) 



Equation (Al), which defines the contribution to the even 
part VPe of the VP at a spatial point (a;', y' , z'), is equivalent 
to the fundamental integral equation for two-integral DFs of 
axisymmetric disc galaxies (e.g. eq. [CI] of HQ). This is not 
unexpected since in both cases the integrations are carried 
out over a two-dimensional velocity space. Equation (C2) 
of HQ transforms the surface "density" of a disc galaxy to 
a pseudo volume-density of an axisymmetric galaxy. Here 
we need its inversion to relate a to the "density" p. It is 



1 [ h 
7r^/2 J V 



($^ 
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in which the subscript 1 denotes the partial derivative with 
respect to the first argument. Completing the z'-integration 
yields 
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VVJvzr,x',y') = Idz' [ ^i|li£2d$, (A4) 
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always provided x'vz' sin i = 0. We need to have the func- 
tion p(^,B?^ to use formula (A4) since it is generally diffi- 
cult to relate /5($,s^) explicitly to the actual density. When 
this function is available explicitly, one simply integrates 
equation (A4). When this function is only know implicitly, 
one must first evaluate /5($,s^) by numerical means and 
then perform numerical integration. 

Below we summarize the three cases in which equation 
(A4) can be used. 



Al Minor axis (a;' = 0) 

Now equation (A4) gives the VP on the projected minor axis 
as 

*(0,y',^')-»=,/2 

VVJvzr,Q,y') = ^ / dz' I 

Pi {—y' cos i + z' sin i)^^ 



^<S{0,y',z')-vl,/2-i 



The minor axis VP is even, hence the mean line-of-sight 
velocity there is identically zero, as is physically evident. 

A2 Face-on (i = 0) 

Now equation (A4) gives the VP of a face-on projection, 
which is completely even, and is given by 

^iR'^.z'^)-^l,/2 

VPe(..,;.',2/')=^j,(^ ^^1^-' J 

^ *i (A6) 
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where K = x + y 



A3 Systemic velocity (■u^/ = 0) 

In this case equation (A4) gives the density of stars with 
vanishing line-of-sight velocity in the whole plane of the 
sky. It is independent of the mean streaming of stars, and 
is given by 

oo ^{x',y',z') 



*o= (A7) 
^<i!{x',y',z')-i 



APPENDIX B: AUXILIARY RESULTS FOR 
SCALE-FREE SPHEROIDS 

The function Pa.,q{d) that appears in the separated form 
(3.11) of the potential of the scale-free spheroids is given by 



oo 

I , \ I f du fsin^e cos= 61 V+^/^l 



when a —2, and by its limit 
du 



-2,9 



1 / du 
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ln(sin=e+4^cos=e), (B2) 
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when a = —2. Here we have defined 
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so that J-2,q = (2/e) arcsin e, with = 1 — <^ 
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To evaluate the function p that appears in equation 
(3.14), we use the explicit formula (3.6) for the density to 
solve for in terms of B? and p, and then substitute the 
result into the expressions (3.7) and (3.8) for the potential. 
This gives 



(a = -2), 
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+2 
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(«/-2), 



(B4) 



where Ja,q is defined above, and 
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For given a and q, expression (B4) can be solved numerically 
to find p for each value of C,, including complex ones. The 
value p{Q) can be found explicitly: 



^(0) 



■ exp{-Kg/I-2,q), a = -2, 



(Jc,,//c,,)"/("+=', a/ -2, 
where we have written 
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so that 1-2, q = J-2,q- lu principle we can also compute 
the derivatives p^^\o) by successive differentiation of the 
integral equations (B4), and so construct a Taylor series for 
p. This process provides an approximation of the function 
p(^^,B?), which consists of elementary density components 
whose DFs can be easily written down (Fricke 1952; Toomre 
1982; HQ; Evans 1994). It is hard in practice to carry out 
this process beyond the first or the second derivative since 
the differentiations soon lead to lengthy expressions, and its 
accuracy is as yet unknown. For instance, a linear approxi- 
mation of /o(^) gives rise to an approximation of /e(7;^) which 
is linear in rf^ , and hence resembles the /e(7;^) of Evans's 
scale-free power-law models. However, for the scale-free 
spheroids this approximation is not accurate unless they are 
very nearly spherical. 

We can compute /e(7;^) accurately by the contour inte- 
gral (2.2). By substitution of equations (3.14) into equation 
(2.2), carrying out partial derivatives, use of equation (3.15), 
transformation to the scaled variables (3.16), and use of cer- 
tain simple transformations, we arrive at integrals for /e(7;^). 
We need to distinguish three cases: a = —2, —3 < a < —2 
and — 2 < a < 0. The result for the a = —2 case is 
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vhere X = rf^ exp(i — l)/i, and the function H is defined 



H{X) = -|^(X) + (2 - |)X^'(X) + X'p"{X). (B9) 



For — 3 < a < — 2 we have 

('0 + ) 



M^')= ^'°~,X'^ J ijf ^,HiY), (BIO) 



where to = 2/(a + 4), and Y = 7;^(i/<o)'°''''°"^'(*o " !)/(*- 
1). When a > -2, 
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where Z = r/^ {-to / tf" - io)/(l + t). When r/ = 
0, X = Y = Z = 0. The reduced contour integrals in 
equations (B8), (BIO) and (Bll) can then be evaluated by 
wrapping the contours tightly around respective branch cuts 
to convert them to elementary real integrals. The results are 



/e(0) 



and 



exp(l)' 



/e(0)=(-")^° >-Wx 
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a < -2, 



2-=^^ 1 5(1,^-1), a>-2, 



(B12) 
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and /o(0) is given in equation (B6). In the numerical evalua- 
tion of integrals (B8) and (Bll) we have employed the new 
integration variables s = exp(i — 1) and s = —1/t, respec- 
tively, so the resulting integrals are along finite paths. 

The case of a density po{m/b')'^ in a non-self-consistent 
power-law potential of the form (3.25), discussed in Sec- 
tion 3.4, leads to DFs of the form (3.28). The functions 
fani^^V'^) that appear there can be computed by integrals 
of the form (B8), (BIO) and (Bll), for cases with 7 = 0, 
7 < and 7 > 0, respectively, with the following modifica- 
tions. We replace a + 2 by 7 in the definition of to. The 
function H now is 



HiX) = ^[C--l)elX + ^-l]il-elxr. 



(B14) 



The exponential term in the integral (B8) for the case 7 = 
is replaced by exp[— a(i — l)/2]. The power terms in the 
integrals (BIO) and (Bll) for the cases 7 have exponent 
2 — (a/7). However, the major difference is that the function 
p{C) is now equal to (1 — e^^)"''^, and hence is explicit, 
so that the integrations can be evaluated numerically in a 
straightforward manner. 



APPENDIX C: VELOCITY DISPERSIONS IN 
THE SCALE-FREE SPHEROIDS 

The explicit solution of the Jeans equations for an f[E,Lz) 
axisymmetric model is (Hunter 1977) 
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Figure Cl. Dynamical quantities in the meridional plane for 
f{E,Lz) singular isothermal spheroids with four different axis 
ratios, as functions of the latitudinal angle 6 (in degrees). The 
value 6 = 0° corresponds to the symmetry axis, the value 6 = 90° 
to the equatorial plane. Solid curves: aji = az] Dashed curves: 
{v^y^^. The quantities are given in units of Vq , the circular 
velocity in the equatorial plane. This figure was adapted from 
van der Marel (1994e). 

The mean streaming motions (vr) and (vz) vanish in an 
axisymmetric model, so that the second moments (v^) and 
(v^) are equal to the dispersions and a^. 

For the scale-free spheroids, the function p(^i,B?) is 
given in terms of the single-variable function p. A simple 
relation between the two identical dispersions ir^, and the 
second moment can then be established by a judicious 
combination of the equations above such that the resulting 
integrand is a total derivative. We find 

{vD - {2a + 3)al, = - (a + 2)*, (C2) 

where $ is the potential given by equation (3.8) (for a ^ —2) 
or (3.7) (for a = —2) and the constant Vo is given in equation 
(3.10). This relation is valid locally, i.e., at each point {R, z). 
Similar relations can be derived for the higher order even 
moments. We note that we must have a < — 1 for (v^) and 
iTjj to be finite everywhere. 

We now restrict our attention to the singular isothermal 
spheroids (i.e., a = —2). The velocity dispersions are then 



elementary. The equation (Cl) for 
as 



IT. can be written 
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Differentiating equation (3.7) gives 
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Upon substitution of 
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the integral in equation (C3) can be carried out, with result: 
1 e ^ ' 

(C6) 
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The remaining dispersion now follows trivially from 

equation (C2), which for a = —2 reduces to 



(C7) 



The constant Vo is equal to the circular velocity in the equa- 
torial plane for these models. Both ir^ and are indepen- 
dent of radius, but do depend on the polar angle 6 defined 
hy R = rsinO and z = rcosO. The total second moment 
parallel to the equatorial plane, (v^) + ir^, is independent 
of 9, and is always equal to Vo . On the z-axis [9 = 0) we 



have aR 



(v^) = Vo /2. The second moments are 



non-negative when < g < 3.46717. This does not imply 
that /e > for all these models. The analysis in Section 
3.2 shows that /e > only when q < 1.3903. Figure Cl 
shows the dynamical quantities of the singular isothermal 
spheroids for various axis ratios, in units of Vo. In the equa- 
torial plane (v^) > a^ for oblate models, and < a^ for 
prolate models. 



APPENDIX D: SEEING CONVOLUTION 

Any observed quantity is a line-of-sight projected quan- 
tity, averaged over some finite pixels on a detector. For 
ground-based observations we also have to take into ac- 
count the effect of atmospheric seeing, which convolves the 
projected properties of the galaxy with a point spread func- 
tion. This PSF is often taken as Gaussian or as the sum of 
Gaussians. Here we restrict ourselves to Gaussians that are 
circular on the plane of the sky. 

Consider a rectangular pixel R of size 21 x 2w, whose 
centre is at the point {x'o,y'o) and whose axes make an ar- 
bitrary angle ©o with the a; '-axis on the sky. We define a 
new (£, jf)-coordinate system with its centre at {x'o,y'o) and 
with axes parallel to the sides of the pixel (Figure Dl). The 
associated transformation is then 

x' = x'o + X cos ©0 , - „, 

(Dl) 



y sin ©0 , 

y' =y'o +x sin Qo +y cos ©0 . 



Let S[x' , y') be a function that depends linearly on the lumi- 
nosity density of the galaxy, e.g., S(a;',2/') x YP[vzi', x' ,y'), 
as defined by equation (2.13). Assuming the PSF to be 
a single normalized Gaussian, Aexp[ — R'^/2a^), the seeing 
convolved function 5s satisfies: 



S,{xo + X ,yo + y ) = A / / dx dy 



oo oo 

II 



(D2) 



The observable So at the position (2:0,2/0) is: 

I w 

So{x'o,y'o) = -^ / S,{x',y')dx'dy' = / dxdy 

JrJ J (D3) 

— I —W 

Ss{xq + X COS ©0 — y sin Oo, 2/0 + ^ sin ©0 -\- y cos 0o). 
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Figure Dl. The coordinate system used for tlie seeing convolu- 
tion of observations with finite pixels of size 2i X 2w. 



Substituting equation (D2) into equation (D3) and exchang- 
ing the order of integrations, we find that the x- and y- 
integrations can be carried out to give expressions in terms 
of error functions. A final use of polar coordinates defined 
hy X = r cos(5 + ©o) and y = r sm[6 + ©o), which contains 
a rotation in the (a;, 2/)-coordinates, yields 



(D4) 



(D5) 



S.{x'o,yo) = ^ J rdr J de K{r,e) 



S[x'o + r cos{d + ©o), y'o + rs\D.{d + ©o)], 

where 

'^("')=^h'(^)-"'(=^)l 

When the PSF is a sum of Gaussians, the kernel K is then 
also a sum of expressions given by equation (D5). Equation 
(D4) therefore allows us to combine seeing convolution and 
pixel averaging in one step. The error functions in the kernel 
K can be computed using efficient algorithms. We note 
that for the given size of pixel and PSF, K is independent 
of the position (a;o,2/o) on the sky, hence can be computed 
beforehand on a grid in {r, 6) space. 
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